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Numerical simulation of high drag reduction 
in a turbulent channel flow with polymer 

additives 

By Yves Dubief 


1. Introduction 

The addition of small amounts of long chain polymer molecules to wall-bounded flows 
can lead to dramatic drag reduction. Although this phenomenon has been known for 
about fifty years, the action of the polymers and its effect on turbulent structures are 
still unclear. Detailed experiments have characterized two distinct regimes (Warholic 
et al. 1999), which are referred to as low drag reduction (LDR) and high drag reduction 
(HDR). The first regime exhibits similar statistical trends as Newtonian flow: the log- 
law region of the mean velocity profile remains parallel to that of the Newtonian flow 
but its lower bound moves away from the wall and the upward shift of the log-region is 
a function of drag reduction, DR. Although streamwise fluctuations are increased and 
transverse ones are reduced, the shape of the rms velocity profiles is not qualitatively 
modified. At higher drag reductions, of the order of 40-50%, the flow enters the HDR 
regime for which the slope of the log-law is dramatically augmented and the Reynolds 
shear stress is small (Warholic et al. 1999; Ptasinski et al. 2001). The drag reduction 
is eventually bounded by a maximum drag reduction (MDR) (Virk & Mickley 1970) 
which is a function of the Reynolds number. While several experiments report mean 
velocity profiles very close to the empirical profile of Virk & Mickley (1970) for MDR 
conditions, the observations regarding the structure of turbulence can differ significantly. 
For instance, Warholic et al. (1999) measured a near-zero Reynolds shear stress, whereas 
a recent experiment (Ptasinski et al. 2001) shows evidence of non-negligible Reynolds 
stress in their MDR flow. To the knowledge of the authors, only the LDR regime has 
been documented in numerical simulations (Sureshkumar et al. 1997; Dimitropoulos et al. 
1998; Min et al. 2001; Dubief & Lele 2001; Sibilla & Baron 2002). This paper discusses 
the simulation of polymer drag reduced channel flow at HDR using the FENE-P (Finite 
Elastic non-linear extensibility-Peterlin) model which was used for the first LDR simu- 
lation by Sureshkumar et al. (1997). Flow and polymer parameters are close to realistic 
polymer drag reducing conditions. High drag reductions are achieved by using finite dif- 
ferences and a robust time stepping technique. A minimal channel flow is also used as 
a numerical experiment to investigate the effect of the outer region turbulent structures 
on the overall drag at HDR. The drag reducing action of the model is finally studied 
through the structure of energy transfers from the polymers to the velocity components. 
This investigation sheds some light on the details of polymer drag reduction. 


2. Governing equations and numerical method 

The formalism of the constitutive equations for viscoelastic flows typically includes 
the assumption of uniform concentration of the polymer solution, and the momentum 
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equations thus become: 
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where j3 is the ratio of the solvent viscosity v to the total viscosity and effectively controls 
the concentration of polymers. The Reynolds number is defined as Re = Uh/ v, where 
the velocity, U, and length, h, scale are defined in the next section. Note the addition of 
the viscoelastic stress which is later referred to as . The viscoelastic tensor in /, is 
obtained by solving the FENE-P equation, 
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where the conformation tensor, Cij, is the phase average of q t qj, qi being the component 
of the end-to-end vector of each individual polymer which has a maximum dimensionless 
extensibility, L. The Weissenberg number We is the ratio of the largest polymer relax- 
ation time A to the flow time scales, such that We = XU/h. The numerical method is 
essentially that of Min et al. (2001) modified to simulate very elastic and long polymer 
molecules. The successful modification, consisting of a novel time advancement scheme 
for Eq. (2.2), is described and validated in Dubief et al. (2003). In the present paper 
only a brief outline of the method is given. The momentum equations are solved on a 
staggered grid with second-order central finite differences. The divergence of the polymer 
stress (Eq. 2. 1) and the spatial derivatives of are computed using a fourth order com- 
pact scheme and a third order upwind compact scheme, respectively. Time advancement 
of Eqs. (2.1) and (2.2) is performed by the classical semi-implicit second-order Crank- 
Nicolson/third-order Runge-Kutta scheme. In the momentum equation, the Newtonian 
viscous stress are treated implicitly in the wall-normal direction. Eq. (2.2) is solved with 
a new semi-implicit time scheme which guarantees the trace of the c l3 remains upper 
bounded ( Ckk < L 2 ). 

In turbulent flows, Eq. (2.2) proves to be fairly stiff and therefore delicate to solve. 
Originally, Sureshkumar et al. (1997) added a diffusive term to the FENE-P ( D(cij ) = 
dkdkCij) and used a fairly significant diffusivity coefficient k. Min et al. (2001) later 
pointed out that the addition of a diffusivity everywhere in the flow causes the smearing of 
polymer stress gradients and suggested a local approach. The Local Artificial Dissipation 
(LAD) is a second-order numerical error ( D(cij ) = A \d\cij, A k denoting the local grid 
size) which is applied only to nodes for which the positiveness of the conformation tensor 
is not satisfied. Min et al. designed the LAD as an-extra diffusion of the scheme used for 
the advection of c,j. 

As discussed in Dubief et al. (2003), the advection term in Eq. (2.2) creates small 
scales similar to the ones observed for a passive scalar at high-Schmidt number, i.e. a 
few orders of magnitude smaller than the dissipative flow scale, the Kolmogorov scale. 
The LAD allows the polymer field to develop sharper gradients than a global diffusivity 
approach, however the quality of the simulation depends strongly on the resolution of 
the smallest turbulent scales. When the drag is very close to its Newtonian value, direct 
numerical simulations can only resolved up to Kolmogorov which is then the cutoff scale 
for the advection of Cij. At HDR, the turbulence is so reduced that the simulation of 
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Figure 1. Left: Normal stress profiles of the viscoelastic simulation L100W120MC, □ , 
compared to the Brownian dynamic simulation , . 


smaller scales for the conformation tensor becomes possible. The quality of a FENE-P 
simulation is assessed by comparing the polymer stress obtained by Eq. (2.2) and the one 
obtained with a Brownian dynamic simulation (Terrapon et al. 2003). Since the Brownian 
simulation is only a one-way coupling from the flow to polymers, the comparison is 
achieved by using the velocity field computed from the viscoelastic simulation. In the 
case of HDR (Fig. 1) the agreement is within 10% which is considered satisfactory. 

Finally, as in any simulation involving dramatic reduction of drag, particular atten- 
tion to the dimensions of the computational domain must be paid. Here, simulations are 
performed in a channel flow, where periodicity is enforced in the homogeneous direction, 
x and z, and no-slip is prescribed on the walls. Periodicity is a very satisfactory bound- 
ary condition providing that the energy-containing turbulent structures are smaller than 
the dimensions of the computational domain. In the case of polymer drag reduction, 
a coarsening of the streaks is observed (White et al. 2003), suggesting that all turbu- 
lent scales in the near-wall region are likely to grow with increasing drag reduction. As 
demonstrated by Jimenez & Moin (1991) in their numerical experiment of the minimal 
flow unit capable of sustaining turbulence, turbulence vanishes when the computational 
domain is smaller than 1000 by 100 wall units in the streamwise and spanwise directions, 
respectively. The wall unit is defined the ratio of the kinematic viscosity v to the friction 
velocity u T = \J v(dU / dy) wa u • Based on the Newtonian friction velocity, the length and 
width of the computational domain have been varied from 1000 to 6000 and from 300 to 
1200, respectively for a drag reduction of DR = 60%. While velocity correlations in the 
streamwise direction do not go exactly to zero, even for Lf = 6000, the difference in drag 
reduction is less than 1% for a domain 4000 x 1200 compared with 6000 x 1200. In the 
spanwise direction, correlations drop to zero for Lf = 1200, which has led us to choose 
the intermediate domain 4000 x 600 x 1200 or Anh x 2h x Ah in integral units, as shown 
in Table 1. In the course of the study of the channel dimensions, it was noticed that tur- 
bulence never disappears in the smallest domain, in spite of its sub-critical dimension in 
drag reduction mode with Lf ~ 600 and L+ ~ 150 (based on the reduced skin-friction). 
Results obtained with this domain are also discussed as a numerical experiment. 

Simulations are performed in a channel flow at an intermediate Reynolds number, 
Re = 7500 based on the half-width h of the channel and the centerline velocity of the 
initial Poiseuille profile. The bulk Reynolds number is ReM = 5000. Conservation of 
the mass flow is imposed which gives h + = hu T /v = 300. The resolution is Ax + = 9, 
A y + = 0.1 — 5 and A z + = 6, when normalized by the skin friction at DR = 0%. Statistics 
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Name 

Line/symbol 

Dimensions 

L 

We T0 

0 

DR 

L60W84LC 

• 

47 ih x 2 h x 4 h 

60 

84 

0.9 

47% 

L60W84MC 

o 

7r h x 2h x h 

60 

84 

0.9 

67% 

L100W120LC 

■ 

47 rh x 2/i x 4A 

100 

120 

0.9 

60% 

L100W120MC 

□ 

nh x 2h x h 

100 

120 

0.9 

72% 


Table 1. Polymer parameters used for the viscoelastic simulations. The Weissenberg number 
We-ro is normalized the wall-shear stress for the Newtonian simulation (DR = 0%). For clarity, 
data from simulation L60W84MC appear only on the mean velocity profile plot (Fig 2). 




Figure 2. Left: Mean velocity profiles scaled with inner variables. : Newtonian simula- 
tion (DR = 0%); : MDR asymptote, U + = 11.71ni/ + + 17. Experimental data (White 

et al. 2003): , DR = 45%; , DR = 67%. Other symbols and lines are defined in 

Table 1. Right: RMS of velocity fluctuations scaled with inner variables. Newtonian simulation 

(DR = 0%): , u ' + ; , v' + . For the viscoelastic simulations, symbols are defined in 

Table 1, u' + is denoted by symbols only; v ,+ is indicated by symbols connected by . 


are collected over 300 to 400 convection times h/U, starting after the transient period, 
typically 200 h/U. 


3. Results 

In the following, the simulations are referred as to by LxxWyyCD, where xx, yy and CD 
designate respectively the dimensionless maximum extensibility L, yy the Weissenberg 
number normalized the Newtonian wall-shear and CD the computational domain, MC 
for minimal channel, LC for large channel (see Table 1). 

Two observations can readily be made from table 1. First the FENE-P model requires 
fairly high Weissenberg numbers to achieve HDR. with long-chain type polymers, while 
HDR has experimentally been measured for We ranging from unity to 100. The need 
for high We can be attributed to the underestimation of polymer stress by the model 
in the case of low extension, and therefore be blamed on the simplicity of the FENE-P 
model, based on a dumbbell, neglecting internal modes. The second observation arises 
from the comparison between drag reductions obtained with the minimal channel and the 
larger domain. For L = 60 and We T 0 = 84, the difference is the largest, with the minimal 
channel showing 20% more reduction than the large domain. As discussed in the previous 
section, turbulence in the minimal channel flow under condition of drag reduction is 
very likely to be affected by the periodic boundary conditions. In the present case, the 
streamwise turbulent intensity is significantly reduced around the centerline (Fig. 2b), 
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Figure 3. Reynolds shear stress (3a) and polymer stress (3b) normalized by u r and v. Symbols 

are defined in Table 1. In Fig. 3b, represents the contribution the polymer stress would 

have to Eq. (3.1) in the case of MDR defined by the velocity profile U + = 11.71n(j/ + ) — 17 and 
— wv + = 0. 

indicating that large scale turbulent structures which are otherwise present in the large 
domain cannot develop. Deprived of this energy coming from the outer region of the flow, 
the near- wall region produces a weaker shear. 

All the mean velocity profiles domain exhibit a change of slope in the log-law region, but 
only the ones obtained with minimal channel flow approach Virk’s asymptote, shown by 
a solid line in Fig. 2a. For L60W84LC (DR = 47%), the slope of the log region is roughly 
twice that of the Newtonian, which ascertains the fact that this simulation belongs to 
the lower part of the HDR. regime, while L100W120LC is intermediate. In the case of 
the minimal channel, increasing We for this simulation had no effect on the amount of 
drag reduction. The simulation L100W120MC seems thus representative of an asymptotic 
state for the constraints imposed by the boundary condition channel flow. Also plotted in 
Fig. 2a are some experimental data obtained at Stanford by White et al. (2003), showing 
that similar results velocity profiles can be obtained in a turbulent boundary layer with 
non-uniform polymer concentration. 

The turbulent velocity fluctuations (Fig. 2b) behave as observed experimentally (see 
White et al. 2003). The peak of u' shifts away from the wall and its magnitude increases 
slowly compared to the Newtonian flow when normalized by u T . The wall-normal com- 
ponent v' follows an opposite trend; w’ behaves as v' and consequently does not appear 
on the plot for clarity. In drag reduced flow, the maximum of u' + is indeed higher or 
comparable with the DR = 0% case, as found in experiments (Warholic et al. 1999; 
Ptasinski et al. 2001; White et al. 2003). The strong reduction of the transverse fluctu- 
ations suggests that polymers target preferentially the vortices, which have a significant 
contribution to the transverse velocity components. 

The balance of stresses, 

( v + \ dU + 

-m+- \l- —j +0—+(1-0 )t+v = 0, (3.1) 

contains significant information regarding the mechanism of polymer interaction with 
the mean flow. As indicated in Fig. 3a, the Reynolds shear stress reduces with increasing 
drag reduction. For DR = 60% (L100W120LC), — uv + is approximately a third of its 
magnitude in the Newtonian case. Conversely, the polymer stress increases with increas- 
ing drag reduction (Fig. 3b). At DR = 60%, its near- wall contribution to Eq. (3.1) has 
the same magnitude as the Reynolds shear stress. In the case of the minimal channel, 
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Figure 4. Snapshot of vortices identified by the Q-criterion (Dubief & Delcayre 2000) and the 
trace ( Ckk/L 2 x 100). Left L60W84MC; right L60W84LC. Vortices are only shown in the lower 
half of the domain. 

the polymer stress is actually 1.5 higher than the Reynolds shear stress, which outlines 
the mechanism by which turbulence is sustained even in a minimal channel flow. The 
level of Reynolds shear stress and polymer stress are comparable to the ones measured 
by Ptasinski et al. (2001), even for their MDR experiment. Also plotted in Fig. 3b is 
the polymer stress contribution in case of zero-Reynolds shear stress as observed in the 
experiments of Warholic et al. (1999) and MDR, i.e. above y + ~ 20 the mean velocity 
profile is assumed to collapse with Virk’s asymptotic profile (see Fig. 2a). In that case, the 
polymer stress is almost twice the one measured in the minimal channel flow for the same 
drag reduction. The different behavior observed in Ptasinski et al. (2001)’s experiment, 
as well as our simulation, and Warholic’s demonstrates the essential role played by the 
polymer stress in the self-sustaining mechanism which produces the MDR turbulence. 

4. Observations and perspectives on the drag reduction mechanism 

The topology of the drag reduced flows is based on the same structures as Newtonian 
wall-turbulence: streaks and vortices. In Fig. 4, showing snapshots of L60W84MC and 
L60W84LC, the streaks are identified by the contours of the polymer extension Ckk/L 2 
at the wall, where Ckk is the trace of c t j . The regions of high and low stretch denotes high 
and low-speed streaks, respectively. The coherence of the streaks is rather impressive 
compared to Newtonian turbulence (not shown here) . The streamwise persistence of the 
streaks has also been observed in White et al. (2003) ’s experiments, as well as their 
spanwise stability. Vortices are identified using isosurfaces of the second invariant Q of 
the velocity gradient tensor. Positive values of Q = — S 2 , where fly = djUi — diUj 

and Sij = djUi + diUj isolate regions where the rotation rate dominates the strain rate, 
which provides a reliable identification method of vortices (Dubief & Delcayre 2000) . The 
minimal channel flow contains a few, elongated vortices while the large domain has a more 
varied population of near- wall vortices, short and long quasi-streamwise vortices and 
hairpin-type vortices. These differences in the topology illustrate the different statistical 
behavior observed between the two simulation. On the side view, polymer extension 
is shown to burst away from the near-wall region intermittently in the upwash and 
downwash flows generated by vortices. 


Numerical simulation of high drag reduction with polymers 


445 




Figure 5. (a) Correlation velocity-viscoelastic stress in the three directions, 

{Ruifi = Uifi/u'if'i, no summation on i). , streamwise direction; , wall-normal; 

, spanwise. (b) Conditional probabilities based on viscoelastic stress fluctuations satisfy- 
ing |/j| < (max ( fi(y )) , — h < y < +h) (the lower curves show the fraction of samples satisfying 

this condition). Upper curves: , P(uf x >0); , P(vf y < 0); , P(wf z < 0). 

Data obtained from simulation L100W120MC. 


The most interesting aspect of polymer drag reduction is to understand how such 
small concentrations of microscopic molecules can drive turbulence to states like the 
ones depicted in Fig. 4. One simple measure of the exchange of energy between the flow 
and polymers is the correlation velocity-viscoelastic stress, 



(4.1) 


as shown in Fig. 5a (the viscoelastic stress /* is defined in Eq. 2.1). The term Uifi 
enters the transport equation for the kinetic energy and through this term only the flow 
can be perturbed by polymers: when pi is positive, polymers tend to enhance velocity 
fluctuations while pi < 0, polymers dampen velocity fluctuations. 

Only results from simulation L100W120MC are plotted since all simulations show the 
same trends. In the near-wall region, streamwise velocity fluctuations and viscoelastic 
stress fluctuations are positively correlated while, for the components normal to the mean 
flow, they are anti-correlated. Similar results were obtained by De Angelis et al. (2002) 
at LDR. The correlations indicates that the viscoelastic stress applies some opposition 
control on the transverse velocity fluctuations and plays a significant role in the increase 
of streamwise velocity fluctuations. Away from the viscous sublayer, the maximum cor- 
relation is about ±0.3, ruling out the possibility that polymers enhance all streamwise 
velocity fluctuations nor oppose all the transverse ones. 

Details of the viscoelastic stress action can only be captured by using conditional 
statistics. In this approach, we are interested in the large fluctuations of the viscoelas- 
tic stress, ft since energy is exchanged between polymers and turbulence through this 
term. The probabilities P((uf x )(y) > 0, y £ [—h,+h]), P((vf y ) < 0 ,y £ [—h,+h]) and 
P((wf z )(y) < 0,y £ [—h,+h]) are computed over samples for which |/j| is larger than 
the maximum of the RMS /' over the entire channel. Also plotted are the fractions of 
samples satisfying such a condition for each component. The probability that uf x > 0 
peaks at 80% around y + = 10, while the probabilities that vf y > 0 and wf z > 0 ex- 
hibit a plateau in between 70% and 80% for y + > 10. These strong probabilities suggest 
that large fluctuations of viscoelastic stress are almost perfectly correlated to either tur- 
bulence enhancement for the streamwise component in the near wall region or to drag 
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Figure 6. Snapshot of vortices (white isosurfaces) and isosurfaces of large viscoelastic stress 
fluctuations: blue, negative fluctuations; red, positive fluctuations. Fig. (a): streamwise viscoelas- 
tic stress; Fig. (b): wall-normal viscoelastic stress. The flow is from left to right, the bottom 
views show the lower half a portion of simulation L100W120LC, at DR = 60%. 


reduction for the transverse components in most of the channel. However it appears that 
large fluctuations of f x away from the wall are highly anti-correlated with u fluctuations, 
indicating that where the mean shear is weak, f x behaves like f y or f z . The separation 
of the drag reducing and drag enhancing activities between velocity components was 
also put forward by the numerical experiments (Dubief & Lele 2001) which isolated the 
action of polymers in the buffer region and also suggested that the drag reducing activ- 
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ity of the FENE-P equations comes predominantly from the wall-normal and spanwise 
components. 

Fig. 6 illustrates the instantaneous location of large fluctuations of f x and f y . The 
view from the side gives an idea of the preferred location of these fluctuations while the 
view from the top indicates the relation between polymers and vortices. As observed 
in the correlation of u and f x , a major portion of the large fluctuations is confined at 
the wall (Fig. 6a). Their occurrence is somewhat related with the presence of quasi- 
streamwise vortices in the buffer region and even above. The near- wall activity of the 
streamwise viscoelastic stress has just been shown (Fig. 5) to enhance streamwise velocity 
fluctuations however, away from the wall, the regions of large fluctuations are attached 
to the vortices in a similar way as the wall-normal viscoelastic stress (Fig. 6b). Large 
wall-normal viscoelastic stress appear to be produced exclusively around vortices and 
the same occurs for the spanwise fluctuations (not shown here). This observation holds 
regardless of the amount of DR. These visualizations show the strong correlation which 
exists between drag reducing polymer action and vortices. 


5. Conclusion 

A set of four viscoelastic simulations has been investigated, all in the high drag re- 
duction regime, i.e change of slope in the log-law of the mean velocity profile and a 
significant decrease of Reynolds shear stress. Numerical experiment were performed in a 
minimal channel to show that the basic topology of drag reduced flow can be contained 
in a small computational domain. In its present form and for the polymer parameters 
under consideration, it appears that the FENE-P model does not produce enough stress 
to damp large scale structures in the outer region of the large domain. As drag reduction 
increases, the polymer stress is found to contribute as much as, and even in the mini- 
mal channel larger than, the Reynolds shear stress to the stress balance in the near-wall 
region. Therefore the polymer stress is essential for sustaining the weak turbulent state 
of MDR. The injection of energy from the polymers to the flow is confined to the very 
near wall region and it affects the streamwise velocity component only. Lastly, it is shown 
that the drag reducing activity targets exclusively the near-wall vortices, main source of 
turbulent friction drag (Kravchenko et al. 1993) and key ingredient of the autonomous 
regeneration cycle of near-wall turbulence (Jimenez & Pinelli 1999). 
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